function [M_t,BI_t, D, DPol, MargDamagesHour, M_bin, pi_bin, Exp_sub_t, Tot_ut]=simulate(params, A_t, Nbar_i, HH_bin, CollPerc_t, Pol_t,Own_t, Inc_t, ...
    Sunit_state, Skwh_state, Scost_state, P_state, P_sale_state, PIns_state, PIns0_state, ...
     Load, RenProd, PPRegion, PPRegion2, PPInt, PPCap, PPInf,PPRandStacked, ...
 ProdMedian,DamagesCoef1,DamagesCoef2, DamagesCoefPoll1,DamagesCoefPoll2, TractSunlightProfile, run_damages)

  

global state_max   statetaxcredit_t statemaxtaxcredit_t sunit_cap_t 
global tract_state_xxwalk_mat ntract 
global  tl_global Sunit_tg Scost_tg Skwh_tg

Sunit_t= sum(tract_state_xxwalk_mat.*(ones(ntract,1)*Sunit_state'),2);
Scost_t= sum(tract_state_xxwalk_mat.*(ones(ntract,1)*Scost_state'),2);
Skwh_t= sum(tract_state_xxwalk_mat.*(ones(ntract,1)*Skwh_state'),2);

if tl_global==1
    Sunit_t= Sunit_tg;
    Scost_t= Scost_tg;
    Skwh_t= Skwh_tg;
end


P_t=sum(tract_state_xxwalk_mat.*(ones(ntract,1)*P_state'),2);
PSale_t=sum(tract_state_xxwalk_mat.*(ones(ntract,1)*P_sale_state'),2);

PIns_t=sum(tract_state_xxwalk_mat.*(ones(ntract,1)*PIns_state'),2);

PIns0_t=sum(tract_state_xxwalk_mat.*(ones(ntract,1)*PIns0_state'),2);

[V_inst, N_i, sub_bin]=calc_V_inst2(params, A_t, Nbar_i,  CollPerc_t, Pol_t, Own_t, Inc_t, ...
    Sunit_t, Scost_t, Skwh_t, P_t, PSale_t, PIns_t, PIns0_t);

Prob_sim_bin=(exp(V_inst))./(1+exp(V_inst));

Prob_sim_bin(V_inst>10)=1;

pi_bin=Prob_sim_bin;

Exp_ut=log(1+exp(V_inst));


M_bin=HH_bin.*N_i.*Prob_sim_bin;

M_t=sum(M_bin,2);

BI_t=sum(Prob_sim_bin.*HH_bin,2);

if (sum(isnan(M_t),"all")>0)
    stop
end



ResProd=calc_res_prod(A_t, M_t, TractSunlightProfile); %yearly res prod in mwh

PPRand=squeeze(PPRandStacked(:,:,1));

D=0;
DPol=0;
MargDamagesHour=0;
if run_damages==1


    [D, DPol, MargDamagesHour]=calc_damages5(ResProd, Load, RenProd, PPRegion, PPRegion2, PPInt, PPCap, PPInf,PPRandStacked, ...
        ProdMedian,DamagesCoef1,DamagesCoef2, DamagesCoefPoll1,DamagesCoefPoll2); % damages\\

end

Exp_sub_bin=pi_bin.*sub_bin;
Exp_sub_t=sum(HH_bin.*Exp_sub_bin,2);

[ sigma, gamma_0, gamma_coll, gamma_pol, gamma_own, gamma_inc, gamma_N, gamma_2N]=unpack_params(params);

Tot_ut=sigma.*sum(sum(HH_bin.*Exp_ut));